DLBCL Heatmap and Ternary Plot Analysis: CD4T-Macrophage Interactions

Figure 5 and Supplementary Figure 7 Generation

Introduction

This analysis generates comprehensive visualizations for DLBCL spatial analysis, focusing on:

• Heatmaps: Pathway enrichment analysis for CD4T cells and Macrophages
• Ternary Plots: Three-way spatial correlation analysis between CD4T, Macrophage, and Tumor cells
• Correlation Analysis: Relationship between EBV status, immune dysfunction, and spatial patterns

The analysis produces figures for the main manuscript (Figure 5) and supplementary materials (Supplementary Figure 7).

Setup and Data Loading

Load Required Libraries and Set Paths

# Set global paths
wdpath <- "/bmbl_data/yuzhou/collaborative/Sizun_lab/INDEPTH/SGCC/SGWT_results/DLBCL_GeoMX/GeoMXNew_code_data_publish/Data/"
results_path <- "/bmbl_data/yuzhou/collaborative/Sizun_lab/INDEPTH/SGCC/SGWT_results/DLBCL_GeoMX/GeoMXNew_code_data_publish/Results"
loadingfunction_wdpath <- "/bmbl_data/yuzhou/collaborative/Sizun_lab/INDEPTH/SGCC/SGWT_results/DLBCL_GeoMX/GeoMXNew_code_data_publish/Code/src/"

# Create results directory if it doesn't exist
if (!dir.exists(results_path)) {
  dir.create(results_path, recursive = TRUE)
}

# Load required libraries
library(Ternary)
library(dplyr)
library(Seurat)
library(GSVA)
library(readr)
library(tidyr)
library(tidyverse)
library(deldir)
library(igraph)
library(progress)
library(qs)
library(ggpubr)
library(SingleCellExperiment)
library(ggrepel)
library(gridExtra)
library(pheatmap)
library(circlize)
library(ComplexHeatmap)
library(knitr)
library(DT)

cat("Libraries loaded successfully!\n")

## Libraries loaded successfully!

cat("Results will be saved to:", results_path, "\n")

## Results will be saved to: /bmbl_data/yuzhou/collaborative/Sizun_lab/INDEPTH/SGCC/SGWT_results/DLBCL_GeoMX/GeoMXNew_code_data_publish/Results

Load Custom Functions and Gene Lists

# Source custom functions
source(file.path(loadingfunction_wdpath, "5-Preload_customized_gene_list.R"))
source(file.path(loadingfunction_wdpath, "4-post_DEG_heatmap_visualization_preloadfunction-mergeMacroTumorT.R"))
source(file.path(loadingfunction_wdpath, "2-rank_SGCC_Impulse_preload_function.R"))
source(file.path(loadingfunction_wdpath, "3-WithinCrossDomainTableAndVisualization_preload_function.R"))

cat("Custom functions loaded successfully!\n")

## Custom functions loaded successfully!

Define Analysis Parameters

# Define cell type pairs for analysis
Celltypepairs <- c("CD4T", "Macro")
SGCC_pair <- "Macro_CD4T"

# Assign cell types
CT.1 <- Celltypepairs[[1]]  # CD4T
CT.2 <- Celltypepairs[[2]]  # Macro

cat("Analysis parameters:\n")

## Analysis parameters:

cat("- Cell type 1 (CT.1):", CT.1, "\n")

## - Cell type 1 (CT.1): CD4T

cat("- Cell type 2 (CT.2):", CT.2, "\n")

## - Cell type 2 (CT.2): Macro

cat("- SGCC pair:", SGCC_pair, "\n")

## - SGCC pair: Macro_CD4T

Load Expression Data and Metadata

# Load expression matrices and metadata
CT.1.exp <- list_exp_matrix_metadata[[CT.1]]
CT.2.exp <- list_exp_matrix_metadata[[CT.2]]

# Create SingleCellExperiment objects
spe_obj_CT1 <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(CT.1.exp),
    logcounts = as.matrix(CT.1.exp)
  ),
  colData = list_exp_matrix_metadata$reformedMeta
)

spe_obj_CT2 <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(CT.2.exp),
    logcounts = as.matrix(CT.2.exp)
  ),
  colData = list_exp_matrix_metadata$reformedMeta
)

cat("Expression data loaded:\n")

## Expression data loaded:

cat("- CT1 (", CT.1, ") dimensions:", dim(CT.1.exp), "\n")

## - CT1 ( CD4T ) dimensions: 18722 149

cat("- CT2 (", CT.2, ") dimensions:", dim(CT.2.exp), "\n")

## - CT2 ( Macro ) dimensions: 18722 149

Data Preprocessing

Load SGCC Matrix and Annotation Data

# Load SGCC matrix
SGCC_matrix <- read.csv(file.path(results_path, "60_SGCC_matrix_DLBCL_mergeMacro_mergeTumor_mergeT.csv"), row.names = 1)

# Create submatrix for relevant cell type combinations
recalled_CTs <- combn(c("Macro", "CD4T", "Tumor"), 2, FUN = function(x) paste(x, collapse = "_"))
SGCC_submatrix <- SGCC_matrix[recalled_CTs, ]

cat("SGCC matrix loaded:\n")

## SGCC matrix loaded:

cat("- Full matrix dimensions:", dim(SGCC_matrix), "\n")

## - Full matrix dimensions: 6 30

cat("- Submatrix dimensions:", dim(SGCC_submatrix), "\n")

## - Submatrix dimensions: 3 30

cat("- Cell type combinations:", paste(recalled_CTs, collapse = ", "), "\n")

## - Cell type combinations: Macro_CD4T, Macro_Tumor, CD4T_Tumor

print(head(SGCC_submatrix))

##               DFCI_12.1   DFCI_13.2    DFCI_14.1  DFCI_15.2  DFCI_17.1
## Macro_CD4T  -0.06397598  0.02582195 -0.050237858  0.1000401 0.22439510
## Macro_Tumor -0.12925790 -0.24821949  0.066721881 -0.1987236 0.10983993
## CD4T_Tumor  -0.07167850 -0.30315113  0.004016523 -0.1662047 0.09466403
##               DFCI_18.2    DFCI_19.2    DFCI_2.2  DFCI_22.2   DFCI_23.2
## Macro_CD4T  -0.06296595 -0.055435479 -0.07189339  0.1902091  0.01876928
## Macro_Tumor -0.12674193  0.002625264 -0.26125648 -0.3338448 -0.08006322
## CD4T_Tumor  -0.59578060 -0.044684406 -0.17899099 -0.1336960 -0.09679291
##                DFCI_3.2    DFCI_4.1   DFCI_7.1     DFCI_8.1 Rochester_11
## Macro_CD4T   0.04956142 -0.03092140  0.1285297  0.002045735   -0.2046232
## Macro_Tumor -0.16276185  0.03231743 -0.3943936 -0.207303874   -0.1885134
## CD4T_Tumor   0.00439033  0.15020420 -0.1713256  0.035439807   -0.2027103
##             Rochester_12 Rochester_13 Rochester_14 Rochester_15 Rochester_16
## Macro_CD4T    -0.1067089   -0.1397385   0.01995796  -0.41699825   -0.4500592
## Macro_Tumor   -0.4169700   -0.2081790  -0.64295481  -0.03447285    0.4399144
## CD4T_Tumor     0.1708963   -0.4484026  -0.34438293  -0.07967041   -0.5670744
##             Rochester_17 Rochester_18 Rochester_19 Rochester_21 Rochester_23
## Macro_CD4T    -0.1822344   0.05913465  0.066418978   -0.1049916   0.01891542
## Macro_Tumor   -0.1691245  -0.45340459 -0.003424294    0.2274223  -0.07026215
## CD4T_Tumor    -0.2246444  -0.11853798  0.022796307    0.1844608  -0.51986793
##             Rochester_25 Rochester_4 Rochester_6  Rochester_7 Rochester_9
## Macro_CD4T   -0.01802326   0.1356985  0.05105050  0.245244323  0.02484695
## Macro_Tumor  -0.03739744  -0.2295247 -0.15476801 -0.312085860 -0.14587496
## CD4T_Tumor    0.25714644  -0.8426158 -0.07738904 -0.003873609 -0.11773520

Process Annotation Data

# Load and process annotation data
annotation <- list_exp_matrix_metadata$annotation

# Standardize annotations
annotation <- annotation %>%
  mutate(Annotation = case_when(
    grepl("CD8", Annotation) ~ "CD8T",
    grepl("CD4", Annotation) ~ "CD4T",
    grepl("Tumor", Annotation) ~ "Tumor",
    grepl("M1", Annotation) ~ "Macro",
    grepl("M2", Annotation) ~ "Macro",
    TRUE ~ Annotation
  ))

# Filter annotations
annotation <- annotation %>%
  filter(!Annotation %in% c("Other", "Neutrophil")) %>%
  filter(!grepl("Tonsil", coreName))

cat("Unique core names after filtering:\n")

## Unique core names after filtering:

print(unique(annotation$coreName))

##  [1] "DFCI_4.1"     "DFCI_8.1"     "DFCI_13.2"    "DFCI_15.2"    "DFCI_17.1"   
##  [6] "DFCI_18.2"    "DFCI_19.2"    "DFCI_22.2"    "DFCI_12.1"    "DFCI_14.1"   
## [11] "DFCI_2.2"     "DFCI_3.2"     "DFCI_7.1"     "DFCI_23.2"    "Rochester_4" 
## [16] "Rochester_6"  "Rochester_11" "Rochester_12" "Rochester_14" "Rochester_18"
## [21] "Rochester_19" "Rochester_21" "Rochester_23" "Rochester_13" "Rochester_15"
## [26] "Rochester_17" "Rochester_25" "Rochester_7"  "Rochester_16" "Rochester_9"

Define Selected Cores

# Specify selected cores for analysis
selected_cores <- c(
  "Rochester_4", "Rochester_6", "Rochester_7", "Rochester_9", 
  "Rochester_11", "Rochester_12", "Rochester_13", "Rochester_14", 
  "Rochester_25", "Rochester_15", "Rochester_16", "Rochester_17", 
  "Rochester_18", "Rochester_19", "Rochester_21", "Rochester_23",
  "DFCI_2.2", "DFCI_3.2", "DFCI_4.1", "DFCI_7.1", "DFCI_8.1",
  "DFCI_12.1", "DFCI_13.2", "DFCI_15.2", "DFCI_17.1",
  "DFCI_18.2", "DFCI_19.2", "DFCI_22.2", "DFCI_23.2"
)

# Filter for selected cores
annotation <- annotation %>%
  filter(coreName %in% selected_cores)

cat("Selected cores for analysis:", length(selected_cores), "\n")

## Selected cores for analysis: 29

cat("Cores after filtering:", length(unique(annotation$coreName)), "\n")

## Cores after filtering: 29

Cell Number Summary

# Summarize cell numbers per core
cellnumber_percore <- annotation %>%
  filter(Annotation %in% c("CD4T", "CD8T", "Macro", "Endothelial")) %>%
  group_by(coreName, Annotation) %>%
  summarise(cell_count = n(), .groups = 'drop') %>%
  mutate(ROI_CT = paste0(coreName, "_", Annotation))

cat("Cell number summary:\n")

## Cell number summary:

summary_stats <- cellnumber_percore %>%
  group_by(Annotation) %>%
  summarise(
    total_cells = sum(cell_count),
    mean_per_core = round(mean(cell_count), 2),
    median_per_core = median(cell_count),
    min_per_core = min(cell_count),
    max_per_core = max(cell_count),
    .groups = 'drop'
  )

DT::datatable(summary_stats, 
              caption = "Cell count statistics by annotation",
              options = list(pageLength = 10, scrollX = TRUE))

Sample Filtering

# Filter samples based on cell number differences
filtered_summary <- filter_samples_summary(
  annotation,
  max_diff = 1,
  cell_types = c(Celltypepairs[[1]], Celltypepairs[[2]]),
  tumor_annotation = "Tumor",
  tumor_upper = 1, 
  tumor_lower = 0,
  logfc = 10
)

# Identify samples to drop
drop.samples <- setdiff(spe_obj_CT1$ROI_rename, selected_cores)

cat("Samples to drop:", length(drop.samples), "\n")

## Samples to drop: 1

if(length(drop.samples) > 0) {
  cat("Dropped samples:", paste(head(drop.samples, 10), collapse = ", "), "\n")
  if(length(drop.samples) > 10) cat("... and", length(drop.samples) - 10, "more\n")
}

## Dropped samples: DFCI_14.1

SGCC Analysis and Metadata Preparation

Prepare SGCC Metadata for Both Cell Types

# Prepare metadata for CT1 (CD4T)
prepare_meta_output_CT1 <- prepare_meta_SGCC(
  spe_obj = spe_obj_CT1,
  SGCC_df = SGCC_matrix,
  cell_pair = c(SGCC_pair),
  CT.use = CT.1,
  condition_name = "EBV_Indicator",
  sample.drop = c(drop.samples),
  batch.factor = "Cohort"
)

# Prepare metadata for CT2 (Macro)
prepare_meta_output_CT2 <- prepare_meta_SGCC(
  spe_obj = spe_obj_CT2,
  SGCC_df = SGCC_matrix,
  cell_pair = c(SGCC_pair),
  CT.use = CT.2,
  condition_name = "EBV_Indicator",
  sample.drop = c(drop.samples),
  batch.factor = "Cohort"
)

cat("SGCC metadata prepared:\n")

## SGCC metadata prepared:

cat("- CT1 samples:", ncol(prepare_meta_output_CT1[[2]]), "\n")

## - CT1 samples: 29

cat("- CT2 samples:", ncol(prepare_meta_output_CT2[[2]]), "\n")

## - CT2 samples: 29

Create Enhanced Metadata

# Create enhanced metadata for CT1 (CD4T)
metaoutput_CT1 <- spe_obj_CT1@colData %>%
  as.data.frame() %>%
  select(ROI_CT, virus_loading, virus_loading_InTumor, total_cells, total_tumor_cells,
         Tumor_other_cell, Tumor_BCL2_cell, Tumor_BCL6_cell, LMP1_g_tumor,CD45RA_CD4T,CD45RO_CD4T,
         Tumor_Myc_cell, LMP1_filtered_mean, LMP1_positive_numebr, LMP1_pct,
         m1_pct, m2_pct, m1_odds, m2_odds, PD1_CD4T, PDL1_macro, Tox_CD4T, 
         LAG3_CD4T, HLADR_macro, dysfunction_score_protein) %>%
  right_join(prepare_meta_output_CT1[[1]], by = "ROI_CT") %>%
  mutate(
    SGCC_cat = case_when(
      Time == 1 ~ "SGCC low",
      Time == 2 ~ "SGCC mediate", 
      Time == 3 ~ "SGCC high",
      TRUE ~ NA_character_
    ),
    EBV_status = case_when(
      Condition == "case" ~ "EBV+",
      Condition == "control" ~ "EBV-"
    )
  )

# Create enhanced metadata for CT2 (Macro)
metaoutput_CT2 <- spe_obj_CT2@colData %>%
  as.data.frame() %>%
  select(ROI_CT, virus_loading, virus_loading_InTumor, total_cells, total_tumor_cells,
         Tumor_other_cell, Tumor_BCL2_cell, Tumor_BCL6_cell, LMP1_g_tumor,CD45RA_CD4T,CD45RO_CD4T,
         Tumor_Myc_cell, LMP1_filtered_mean, LMP1_positive_numebr, LMP1_pct,
         m1_pct, m2_pct, m1_odds, m2_odds, PD1_CD4T, PDL1_macro, Tox_CD4T, 
         LAG3_CD4T, HLADR_macro, dysfunction_score_protein) %>%
  right_join(prepare_meta_output_CT2[[1]], by = "ROI_CT") %>%
  mutate(
    SGCC_cat = case_when(
      Time == 1 ~ "SGCC low",
      Time == 2 ~ "SGCC mediate",
      Time == 3 ~ "SGCC high", 
      TRUE ~ NA_character_
    ),
    EBV_status = case_when(
      Condition == "case" ~ "EBV+",
      Condition == "control" ~ "EBV-"
    )
  )

cat("Enhanced metadata created:\n")

## Enhanced metadata created:

cat("- CT1 metadata rows:", nrow(metaoutput_CT1), "\n")

## - CT1 metadata rows: 29

cat("- CT2 metadata rows:", nrow(metaoutput_CT2), "\n")

## - CT2 metadata rows: 29

cat("- EBV+ samples:", sum(metaoutput_CT1$EBV_status == "EBV+", na.rm = TRUE), "\n")

## - EBV+ samples: 17

cat("- EBV- samples:", sum(metaoutput_CT1$EBV_status == "EBV-", na.rm = TRUE), "\n")

## - EBV- samples: 12

GSVA Pathway Analysis

Run GSVA for CD4T Cells

# Set up GSVA parameters for CD4T
SetGSVAPar_CT1 <- gsvaParam(
  exprData = prepare_meta_output_CT1[[2]],
  geneSets = functionalenrichment[grep(paste0("^", Celltypepairs[[1]]), names(functionalenrichment))],
  minSize = 3, 
  maxSize = 400,
  kcdf = "Gaussian"
)

# Run GSVA analysis
gsva_results_CT1 <- gsva(SetGSVAPar_CT1, verbose = TRUE)
GSVA_pathwayname_CT1 <- rownames(gsva_results_CT1)

cat("GSVA analysis for CD4T completed:\n")

## GSVA analysis for CD4T completed:

cat("- Number of pathways:", nrow(gsva_results_CT1), "\n")

## - Number of pathways: 12

cat("- Pathways analyzed:", paste(GSVA_pathwayname_CT1, collapse = ", "), "\n")

## - Pathways analyzed: CD4T_T_Cell_Proliferation, CD4T_Antigen_Receptor_signaling, CD4T_Regulatory_T_Cell_Differentiation, CD4T_Programmed_Cell_Death, CD4T_Cellular_Senescence, CD4T_Cell_Population_Proliferation, CD4T_T_Cell_Mediated_Immunity, CD4T_Regulation_Of_Viral_Process, CD4T_Positive_Regulation_Of_TNF, CD4T_PancancerT_Adhesion, CD4T_customized_Tcell_exhaustion, CD4T_Innate_Immune_System

Run GSVA for Macrophages

# Set up GSVA parameters for Macrophages
SetGSVAPar_CT2 <- gsvaParam(
  exprData = prepare_meta_output_CT2[[2]],
  geneSets = functionalenrichment[grep(paste0("^", Celltypepairs[[2]]), names(functionalenrichment))],
  minSize = 3, 
  maxSize = 100,
  kcdf = "Gaussian"
)

# Run GSVA analysis
gsva_results_CT2 <- gsva(SetGSVAPar_CT2, verbose = TRUE)
GSVA_pathwayname_CT2 <- rownames(gsva_results_CT2)

cat("GSVA analysis for Macrophages completed:\n")

## GSVA analysis for Macrophages completed:

cat("- Number of pathways:", nrow(gsva_results_CT2), "\n")

## - Number of pathways: 21

cat("- Pathways analyzed:", paste(GSVA_pathwayname_CT2, collapse = ", "), "\n")

## - Pathways analyzed: Macro_merged_translation, Macro_Macrophage_differentiation, Macro_Gene_expression, Macro_TOR_signaling, Macro_Regulation_T_Cell_Differentiation, Macro_MHC_class_II, Macro_NF_kappaB, Macro_Response_To_Hypoxia, Macro_Negative_Regulation_Of_Inflammatory_Response, Macro_Signal_Transduction_In_Absence_Of_Ligand, Macro_Negative_Regulation_Of_Tumor_Necrosis_Factor_Production, Macro_Positive_Regulation_Of_T_Cell_Activation, Macro_MHC_Class_II_Protein_Complex_Assembly, Macro_customize_M2, Macro_customize_M1, Macro_MoMac_VERSE_cluster_2, Macro_MoMac_VERSE_cluster_3, Macro_MoMac_VERSE_cluster_6, Macro_MoMac_VERSE_cluster_13, Macro_MoMac_VERSE_cluster_16, Macro_MoMac_VERSE_cluster_17

Heatmap Generation

Generate Macrophage Pathway Heatmap

# Create macrophage pathway heatmap
Genelist_CT2 <- intersect(Macro_customized_M1M2, rownames(prepare_meta_output_CT2[[2]]))

Macropage_heatmap_pathway <- create_complex_heatmap_FixEBV(
  gene_list = NA,
  pathway_matrix = gsva_results_CT2,
  scale_color_ht = 2,
  expr_data = prepare_meta_output_CT2[[2]],
  meta_data = metaoutput_CT2,
  cluster_rows = TRUE,
  cluster_columns = FALSE
)

# Display the heatmap
draw(Macropage_heatmap_pathway)

# Save the heatmap
pdf(file.path(results_path, "Figures_5A_5B_SuppFig7A_Macrophage-CD4T-Macrophage_heatmap.pdf"), 
    width = 10, height = 9)
draw(Macropage_heatmap_pathway)
dev.off()

## png 
##   2

cat("Macrophage heatmap saved to:", file.path(results_path, "Figures_5A_5B_SuppFig7A_Macrophage-CD4T-Macrophage_heatmap.pdf"), "\n")

## Macrophage heatmap saved to: /bmbl_data/yuzhou/collaborative/Sizun_lab/INDEPTH/SGCC/SGWT_results/DLBCL_GeoMX/GeoMXNew_code_data_publish/Results/Figures_5A_5B_SuppFig7A_Macrophage-CD4T-Macrophage_heatmap.pdf

Generate CD4T Pathway Heatmap

# CD4 gene heatmap (based on 05_figure5A&B_SuppFig7A&B_half.R)
Genelist_CT1 <- intersect(Tumor_customized_M1M2_formation, rownames(prepare_meta_output_CT1[[2]]))

# CD4 pathway heatmap (based on 05_figure5A&B_SuppFig7A&B_half.R)
CD4T_heatmap_pathway <- create_complex_heatmap_FixEBV(
  gene_list = NA,
  scale_color_ht = 2,
  pathway_matrix = gsva_results_CT1,
  expr_data = prepare_meta_output_CT1[[2]],
  meta_data = metaoutput_CT1,
  cluster_rows = TRUE,
  cluster_columns = FALSE
)

# Display the heatmap
draw(CD4T_heatmap_pathway)

# Save Tumor pathway heatmap
pdf(file.path(results_path, "Figures_5A_5B_SuppFig7A_CD4T-CD4T-Macrophage_pathway_heatmap.pdf"), 
    width = 10, height = 9)
draw(CD4T_heatmap_pathway)
dev.off()

## png 
##   2

cat("Tumor pathway heatmap saved to:", file.path(results_path, "Figures_5A_5B_SuppFig7A_CD4T-CD4T-Macrophage_pathway_heatmap.pdf"), "\n")

## Tumor pathway heatmap saved to: /bmbl_data/yuzhou/collaborative/Sizun_lab/INDEPTH/SGCC/SGWT_results/DLBCL_GeoMX/GeoMXNew_code_data_publish/Results/Figures_5A_5B_SuppFig7A_CD4T-CD4T-Macrophage_pathway_heatmap.pdf

Generate CD4T Exhaustion Score Annotation

# Create CD4T exhaustion score annotation heatmap
CD4T_heatmap_annotation_exhaustion_score <- create_complex_heatmap_annotation_from_pathway(
  pathway_matrix = gsva_results_CT1,
  scale_color_ht = 2,
  meta_data = metaoutput_CT1,
  pathwayname = "CD4T_customized_Tcell_exhaustion"
)

# Display the annotation heatmap
draw(CD4T_heatmap_annotation_exhaustion_score)

# Save the annotation heatmap
pdf(file.path(results_path, "Figure5B_T_dysfunction_RNA_annotation_bar.pdf"), 
    width = 10, height = 9)
draw(CD4T_heatmap_annotation_exhaustion_score)
dev.off()

## png 
##   2

cat("CD4T exhaustion annotation saved to:", file.path(results_path, "Figure5B_T_dysfunction_RNA_annotation_bar.pdf"), "\n")

## CD4T exhaustion annotation saved to: /bmbl_data/yuzhou/collaborative/Sizun_lab/INDEPTH/SGCC/SGWT_results/DLBCL_GeoMX/GeoMXNew_code_data_publish/Results/Figure5B_T_dysfunction_RNA_annotation_bar.pdf

Ternary Plot Analysis

Prepare Ternary Data

minmax_rescale <- function(x) (x - min(x)) / (max(x) - min(x))
transformed <- apply(as.matrix(SGCC_submatrix), 1, FUN = minmax_rescale)
transformed <- transformed / rowSums(transformed) * 100
transformed <- cbind.data.frame(transformed, t(SGCC_submatrix))
colnames(transformed) <- c("T_M", "T_Tu", "M_Tu", "SGCC_T_M", "SGCC_T_Tu", "SGCC_M_Tu")
transformed$Core <- rownames(transformed)
transformed$ROI_CT <- paste0(rownames(transformed), "_CD4T")

# Join EBV_status (and other metadata) onto ternary data
transformed <- transformed %>% right_join(metaoutput_CT1, by = c("ROI_CT" = "ROI_CT"))
transformed <- transformed %>% filter(!is.na(T_M))

Define Ternary Plot Helper Functions

# Create basic ternary plot using Ternary package
create_ternary_plot_basic <- function(data, title = "Ternary Plot", show_density = FALSE, density_resolution = 30L) {
  # Set up colors for EBV status
  ebv_colors <- c("EBV+" = "#ca6938", "EBV-" = "#4a9d7a")
  ebv_shapes <- c("EBV+" = 21, "EBV-" = 24)
  
  # Create the ternary plot
  TernaryPlot(alab = "M-Tu", blab = "T-Tu", clab = "T-M", 
              main = title, grid.minor.lines = 0)

  # Optionally overlay EBV-specific density contours
  if (show_density) {
    for (status in c("EBV+", "EBV-")) {
      subset_data <- data[data$EBV_status == status, c("M_Tu", "T_Tu", "T_M")]
      if (nrow(subset_data) > 2 && all(colSums(is.na(subset_data)) < nrow(subset_data))) {
        # Remove rows with any NA to avoid errors in density computation
        coords <- as.matrix(subset_data[stats::complete.cases(subset_data), , drop = FALSE])
        if (nrow(coords) > 2) {
          TernaryDensityContour(coords, resolution = density_resolution, color.palette = ebv_colors[status], lwd = 1.5)
        }
      }
    }
  }
  
  # Add points for each EBV status
  for (status in c("EBV+", "EBV-")) {
    subset_data <- data[data$EBV_status == status, ]
    if (nrow(subset_data) > 0) {
      TernaryPoints(subset_data[, c("M_Tu", "T_Tu", "T_M")], 
                   col = ebv_colors[status], 
                   pch = ebv_shapes[status], 
                   cex = 2.5, 
                   bg = ebv_colors[status])
    }
  }
  
  # Add legend
  legend("topright", legend = names(ebv_colors), 
         col = ebv_colors, pch = ebv_shapes, pt.bg = ebv_colors,
         title = "EBV Status", bty = "n")
}

# Create ternary plot with labels
create_ternary_plot_with_labels <- function(data, title = "Ternary Plot with Labels", show_density = FALSE, density_resolution = 30L) {
  # Set up colors for EBV status
  ebv_colors <- c("EBV+" = "#ca6938", "EBV-" = "#4a9d7a")
  ebv_shapes <- c("EBV+" = 21, "EBV-" = 24)
  
  # Create the ternary plot
  TernaryPlot(alab = "M-Tu", blab = "T-Tu", clab = "T-M", 
              main = title, grid.minor.lines = 0)

  # Optionally overlay EBV-specific density contours
  if (show_density) {
    for (status in c("EBV+", "EBV-")) {
      subset_data <- data[data$EBV_status == status, c("M_Tu", "T_Tu", "T_M")]
      if (nrow(subset_data) > 2 && all(colSums(is.na(subset_data)) < nrow(subset_data))) {
        coords <- as.matrix(subset_data[stats::complete.cases(subset_data), , drop = FALSE])
        if (nrow(coords) > 2) {
          TernaryDensityContour(coords, resolution = density_resolution, color.palette = ebv_colors[status], lwd = 1.5)
        }
      }
    }
  }
  
  # Add points for each EBV status
  for (status in c("EBV+", "EBV-")) {
    subset_data <- data[data$EBV_status == status, ]
    if (nrow(subset_data) > 0) {
      TernaryPoints(subset_data[, c("M_Tu", "T_Tu", "T_M")], 
                   col = ebv_colors[status], 
                   pch = ebv_shapes[status], 
                   cex = 2.5, 
                   bg = ebv_colors[status])
      
      # Add text labels
      TernaryText(subset_data[, c("M_Tu", "T_Tu", "T_M")], 
                 labels = subset_data$ROI_CT, 
                 cex = 0.4, 
                 col = "black")
    }
  }
  
  # Add legend
  legend("topright", legend = names(ebv_colors), 
         col = ebv_colors, pch = ebv_shapes, pt.bg = ebv_colors,
         title = "EBV Status", bty = "n")
}

# Create function for biomarker ternary plots
create_biomarker_ternary_plot <- function(data, biomarker_col, title, 
                                         low_color = "#fafafa", high_color = "#fc7051",
                                         show_density = FALSE, density_resolution = 30L) {
  # Set up colors and shapes for EBV status
  ebv_colors <- c("EBV+" = "#ca6938", "EBV-" = "#4a9d7a")
  ebv_shapes <- c("EBV+" = 21, "EBV-" = 24)
  
  # Create the ternary plot
  TernaryPlot(alab = "M-Tu", blab = "T-Tu", clab = "T-M", 
              main = title, grid.minor.lines = 0)
  
  # Optionally overlay EBV-specific density contours
  if (show_density) {
    for (status in c("EBV+", "EBV-")) {
      subset_data <- data[data$EBV_status == status, c("M_Tu", "T_Tu", "T_M")]
      if (nrow(subset_data) > 2 && all(colSums(is.na(subset_data)) < nrow(subset_data))) {
        coords <- as.matrix(subset_data[stats::complete.cases(subset_data), , drop = FALSE])
        if (nrow(coords) > 2) {
          TernaryDensityContour(coords, resolution = density_resolution, col = ebv_colors[status], lwd = 1.5)
        }
      }
    }
  }
  
  # Normalize biomarker values for color mapping
  biomarker_values <- data[[biomarker_col]]
  biomarker_norm <- (biomarker_values - min(biomarker_values, na.rm = TRUE)) / 
                   (max(biomarker_values, na.rm = TRUE) - min(biomarker_values, na.rm = TRUE))
  
  # Create color gradient
  biomarker_colors <- colorRampPalette(c(low_color, high_color))(100)
  point_colors <- biomarker_colors[pmax(1, pmin(100, round(biomarker_norm * 99 + 1)))]
  
  # Add points for each EBV status
  for (status in c("EBV+", "EBV-")) {
    subset_indices <- which(data$EBV_status == status)
    if (length(subset_indices) > 0) {
      subset_data <- data[subset_indices, ]
      TernaryPoints(subset_data[, c("M_Tu", "T_Tu", "T_M")], 
                   col = "black", 
                   pch = ebv_shapes[status], 
                   cex = 2.5, 
                   bg = point_colors[subset_indices])
    }
  }
  
  # Add legend for EBV status
  legend("topright", legend = names(ebv_shapes), 
         col = "black", pch = ebv_shapes, pt.bg = ebv_colors,
         title = "EBV Status", bty = "n")
  
  # Add color bar legend for biomarker (simplified)
  legend("topleft", legend = c("Low", "High"), 
         col = "black", pch = 21, pt.bg = c(low_color, high_color),
         title = biomarker_col, bty = "n")
}

# Create function for AES ternary plots (with log transformation and different colors)
create_aes_ternary_plot <- function(data, aes_col, title, show_density = FALSE, density_resolution = 30L) {
  # Transform the AES data
  data_copy <- data
  data_copy[[paste0(aes_col, "_log")]] <- log1p(data[[aes_col]])
  
  # Create the plot with blue-yellow color scheme
  create_biomarker_ternary_plot(data_copy, paste0(aes_col, "_log"), title, 
                               low_color = "#2424db", high_color = "yellow",
                               show_density = show_density, density_resolution = density_resolution)
}

Basic Ternary Plots

# Basic ternary plot using Ternary package
create_ternary_plot_basic(transformed, "CD4T-Macrophage-Tumor Spatial Correlations")

# Ternary plot with sample names
create_ternary_plot_with_labels(transformed, "Ternary Plot with Sample Labels")

Protein Marker Ternary Plots

# LAG3 CD4T ternary plot
create_biomarker_ternary_plot(transformed, "LAG3_CD4T", "LAG3 Expression in CD4T Cells")

# TOX CD4T ternary plot
create_biomarker_ternary_plot(transformed, "Tox_CD4T", "TOX Expression in CD4T Cells")

# HLA-DR Macrophage ternary plot
create_biomarker_ternary_plot(transformed, "HLADR_macro", "HLA-DR Expression in Macrophages")

# PD-L1 Macrophage ternary plot
create_biomarker_ternary_plot(transformed, "PDL1_macro", "PD-L1 Expression in Macrophages")

Load AES Data and Create Enhanced Ternary Plots

# Load AES (Adaptive Enrichment Score) results
CD4T_Macro <- qread(file.path(wdpath, "SGCC_relevant_analysis", "Macro_CD4Tstat.result.qs"))
Tumor_Macro <- qread(file.path(wdpath, "SGCC_relevant_analysis", "Macro_Tumorstat.result.qs"))
CD4T_Tumor <- qread(file.path(wdpath, "SGCC_relevant_analysis", "CD4T_Tumorstat.result.qs"))

# Rename columns and add Core identifier
colnames(CD4T_Tumor) <- paste0(colnames(CD4T_Tumor), "_T_Tu")
CD4T_Tumor$Core <- rownames(CD4T_Tumor)

colnames(Tumor_Macro) <- paste0(colnames(Tumor_Macro), "_M_Tu")
Tumor_Macro$Core <- rownames(Tumor_Macro)

colnames(CD4T_Macro) <- paste0(colnames(CD4T_Macro), "_T_M")
CD4T_Macro$Core <- rownames(CD4T_Macro)

# Add exhaustion score
sorted_CD4T_customized_Tcell_exhaustion <- gsva_results_CT1["CD4T_customized_Tcell_exhaustion", ][transformed$Sample]

# Create comprehensive dataset
transformed_CD4T_Macro_Tumor <- Tumor_Macro %>%
  filter(Core %in% selected_cores) %>%
  left_join(CD4T_Macro, by = "Core") %>%
  left_join(CD4T_Tumor, by = "Core") %>%
  left_join(transformed, by = "Core")

transformed_CD4T_Macro_Tumor$CD4T_exhaustion_gene <- sorted_CD4T_customized_Tcell_exhaustion

cat("AES data loaded and merged:\n")

## AES data loaded and merged:

cat("- Final dataset dimensions:", dim(transformed_CD4T_Macro_Tumor), "\n")

## - Final dataset dimensions: 29 100

EBV+ vs EBV- Statistical Tests for Scores

# scores to test (only keep existing columns)
score_cols <- intersect(c(
  "LAG3_CD4T","Tox_CD4T","HLADR_macro","PDL1_macro","PD1_CD4T","CD45RO_CD4T",
  "LMP1_filtered_mean","LMP1_pct",
  "CD4T_exhaustion_gene","dysfunction_score_protein",
  "AES_M_Tu","AES_T_M","AES_T_Tu",
  "m2_pct"
), colnames(transformed_CD4T_Macro_Tumor))

ebv_test_one <- function(df, sc) {
  d <- df %>%
    transmute(
      EBV_status = factor(EBV_status, levels = c("EBV-","EBV+")),
      Cohort = factor(Batch),
      value = .data[[sc]]
    ) %>%
    filter(!is.na(EBV_status), !is.na(value), !is.na(Cohort))

  fit <- lm(value ~ EBV_status + Cohort, data = d)
  tibble(
    score = sc,
    estimate = unname(coef(fit)[["EBV_statusEBV+"]]),
    p_value  = summary(fit)$coefficients["EBV_statusEBV+", "Pr(>|t|)"]
  )
}

ebv_score_tests <- bind_rows(lapply(score_cols, \(sc)
  ebv_test_one(transformed_CD4T_Macro_Tumor, sc)
)) %>% select(score, p_value) %>% 
  arrange(p_value)

ebv_score_tests

Enhanced Ternary Plots with AES and Gene Expression

# LMP1 expression ternary plot
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "LMP1_filtered_mean", "LMP1 Expression (Filtered Mean)")

# PD-1 CD4T ternary plot
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "PD1_CD4T", "PD-1 Expression in CD4T Cells")

# CD4T exhaustion gene signature
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "CD4T_exhaustion_gene", "CD4T Exhaustion Gene Signature")

# Dysfunction score protein
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "dysfunction_score_protein", "Dysfunction Score (Protein)")

AES Score Ternary Plots

# AES Macrophage-Tumor
create_aes_ternary_plot(transformed_CD4T_Macro_Tumor, "AES_M_Tu", "AES Score: Macrophage-Tumor")

# AES CD4T-Macrophage
create_aes_ternary_plot(transformed_CD4T_Macro_Tumor, "AES_T_M", "AES Score: CD4T-Macrophage")

# AES CD4T-Tumor
create_aes_ternary_plot(transformed_CD4T_Macro_Tumor, "AES_T_Tu", "AES Score: CD4T-Tumor")

Additional Biomarker Ternary Plots

# LMP1 tumor percentage
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "LMP1_pct", "LMP1 Percentage in Tumor")

# M2 percentage
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "m2_pct", "M2 Macrophage Percentage")

Save Comprehensive Ternary Plot Figure

# Save comprehensive ternary plots using base plotting layout
pdf(file.path(results_path, "Figure_5C_5D_SuppFig7C_Comprehensive_Ternary_Plots.pdf"), 
    width = 30, height = 30)

# 3 rows x 5 columns
par(mfrow = c(3, 5), mar = c(2, 2, 3, 2))

# Row 1
create_ternary_plot_basic(transformed_CD4T_Macro_Tumor, "CD4T-Macrophage-Tumor Spatial Correlations")
create_ternary_plot_with_labels(transformed_CD4T_Macro_Tumor, "Ternary Plot with Sample Labels")
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "CD4T_exhaustion_gene", "CD4T Exhaustion Gene Signature")
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "dysfunction_score_protein", "Dysfunction Score (Protein)")
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "LMP1_filtered_mean", "LMP1 Expression (Filtered Mean)")

# Row 2
create_aes_ternary_plot(transformed_CD4T_Macro_Tumor, "AES_T_M", "AES Score: CD4T-Macrophage")
create_aes_ternary_plot(transformed_CD4T_Macro_Tumor, "AES_M_Tu", "AES Score: Macrophage-Tumor")
create_aes_ternary_plot(transformed_CD4T_Macro_Tumor, "AES_T_Tu", "AES Score: CD4T-Tumor")
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "HLADR_macro", "HLA-DR Expression in Macrophages")
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "PDL1_macro", "PD-L1 Expression in Macrophages")

# Row 3
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "Tox_CD4T", "TOX Expression in CD4T Cells")
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "LAG3_CD4T", "LAG3 Expression in CD4T Cells")
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "PD1_CD4T", "PD-1 Expression in CD4T Cells")
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "LMP1_pct", "LMP1 Percentage in Tumor")
create_biomarker_ternary_plot(transformed_CD4T_Macro_Tumor, "m2_pct", "M2 Macrophage Percentage")

dev.off()

## png 
##   2

# Reset plotting parameters
par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

cat("Comprehensive ternary plots saved to:", file.path(results_path, "Figure_5C_5D_SuppFig7C_Comprehensive_Ternary_Plots.pdf"), "\n")

## Comprehensive ternary plots saved to: /bmbl_data/yuzhou/collaborative/Sizun_lab/INDEPTH/SGCC/SGWT_results/DLBCL_GeoMX/GeoMXNew_code_data_publish/Results/Figure_5C_5D_SuppFig7C_Comprehensive_Ternary_Plots.pdf

write.csv(transformed_CD4T_Macro_Tumor, file =file.path(results_path, "tableSupp4_SGCCternary.csv") )

Correlation Analysis

Prepare Data for Correlation Analysis

# Define functional variables for LMP1 correlation analysis
functional_test_for_LMP1 <- c(
  "dysfunction_score_protein", "LMP1_filtered_mean", "LMP1_g_tumor", 
  "LMP1_positive_numebr", "m2_pct", "m1_pct", "CD4T_exhaustion_gene"
)

# Select data for correlation analysis
selected_data <- transformed_CD4T_Macro_Tumor %>%
  select(all_of(functional_test_for_LMP1))

cat("Correlation analysis variables:\n")

## Correlation analysis variables:

cat(paste(functional_test_for_LMP1, collapse = ", "), "\n")

## dysfunction_score_protein, LMP1_filtered_mean, LMP1_g_tumor, LMP1_positive_numebr, m2_pct, m1_pct, CD4T_exhaustion_gene

cat("Data dimensions for correlation:", dim(selected_data), "\n")

## Data dimensions for correlation: 29 7

Compute Correlations and Generate Scatter Plots

# Compute Spearman correlation matrix
correlation.tst <- cor(selected_data, use = "pairwise.complete.obs", method = "spearman")
correlation.tst <- ifelse(is.na(correlation.tst), 0, correlation.tst)
diag(correlation.tst) <- 0

# Initialize p-value matrix
correlation_pvalues <- matrix(NA, ncol = ncol(selected_data), nrow = ncol(selected_data))
rownames(correlation_pvalues) <- colnames(selected_data)
colnames(correlation_pvalues) <- colnames(selected_data)

# Perform pairwise correlation tests and generate scatter plots
significant_pairs <- 0
for (i in seq_len(ncol(selected_data))) {
  for (j in seq_len(ncol(selected_data))) {
    if (i != j) {
      # Perform Spearman correlation test
      test_result <- cor.test(selected_data[[i]], selected_data[[j]], method = "spearman")
      correlation_pvalues[i, j] <- test_result$p.value
      
      if(correlation_pvalues[i, j] <= 0.05) {
        significant_pairs <- significant_pairs + 1
        
        # Generate scatter plot with ggplot2
        plot_data <- data.frame(
          x = selected_data[[i]], 
          y = selected_data[[j]], 
          EBV = transformed_CD4T_Macro_Tumor$EBV_status
        )
        
        plot <- ggplot(plot_data, aes(x = x, y = y, color = EBV)) +
          geom_point(size = 3) +
          geom_smooth(method = "lm", color = "blue", se = TRUE) +
          scale_color_manual(
            values = c("#ca6938", "#4a9d7a"),
            breaks = c("EBV+", "EBV-")
          ) +
          theme_classic() +
          labs(
            title = paste(colnames(selected_data)[i], "vs", colnames(selected_data)[j]),
            x = colnames(selected_data)[i],
            y = colnames(selected_data)[j]
          ) + # geom_smooth(span = 0.3)+
          theme(
            panel.background = element_rect(fill = "transparent", color = NA),
            plot.background = element_rect(fill = "transparent", color = NA),
            legend.background = element_rect(fill = "transparent", color = NA)
          )
        
        # Save scatter plot
        ggsave(
          filename = file.path(results_path, paste0("SuppFig7_scatterplot_", 
                                                   colnames(selected_data)[i], "_vs_", 
                                                   colnames(selected_data)[j], ".pdf")),
          plot = plot,
          device = "pdf",
          bg = "transparent",
          width = 8, height = 8, units = "in", dpi = 300
        )
      }
    }
  }
}

cat("Correlation analysis completed:\n")

## Correlation analysis completed:

cat("- Significant pairs (p ≤ 0.05):", significant_pairs, "\n")

## - Significant pairs (p ≤ 0.05): 16

cat("- Scatter plots saved to Results directory\n")

## - Scatter plots saved to Results directory

Generate Correlation Heatmap

# Create indicator matrix for significance
indicator_matrix <- ifelse(correlation_pvalues < 0.05 | is.na(correlation_pvalues), 1, 0.3)
heatmap_mat_indicator <- correlation.tst * indicator_matrix

cols <- colorRampPalette(c("#528ad0", "white", "#cb438d"))(101)

alph_non_signaficant <- pheatmap(
  heatmap_mat_indicator,
  color = cols,
  breaks = seq(-1, 1, length.out = length(cols) + 1),
  na_col = "#ffffff",
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  show_rownames = TRUE,
  show_colnames = TRUE,
  main = "Correlation Matrix (Significant Correlations Highlighted)"
)


# Display the heatmap
print(alph_non_signaficant)

Correlation Results Summary

# Create summary of significant correlations
sig_correlations <- data.frame()
for (i in seq_len(nrow(correlation_pvalues))) {
  for (j in seq_len(ncol(correlation_pvalues))) {
    if (i != j && !is.na(correlation_pvalues[i, j]) && correlation_pvalues[i, j] <= 0.05) {
      sig_correlations <- rbind(sig_correlations, data.frame(
        Variable1 = rownames(correlation_pvalues)[i],
        Variable2 = colnames(correlation_pvalues)[j],
        Correlation = round(correlation.tst[i, j], 4),
        P_value = round(correlation_pvalues[i, j], 4)
      ))
    }
  }
}

if(nrow(sig_correlations) > 0) {
  cat("Significant correlations found:\n")
  DT::datatable(sig_correlations, 
                caption = "Significant correlations (p ≤ 0.05)",
                options = list(pageLength = 10, scrollX = TRUE)) %>%
    DT::formatRound(columns = c("Correlation", "P_value"), digits = 4)
} else {
  cat("No significant correlations found at p ≤ 0.05\n")
}

## Significant correlations found:

Files Generated

cat("=== OUTPUT FILES (ternary / correlation) ===\n")

## === OUTPUT FILES (ternary / correlation) ===

output_files <- c(
  "Figure_5C_5D_SuppFig7C_Comprehensive_Ternary_Plots.pdf",
  "SuppFig7_LMP1_correlation_heatmap.pdf"
)

for (f in output_files) {
  fp <- file.path(results_path, f)
  cat(if (file.exists(fp)) "✓" else "✗", f, "\n")
}

## ✓ Figure_5C_5D_SuppFig7C_Comprehensive_Ternary_Plots.pdf 
## ✗ SuppFig7_LMP1_correlation_heatmap.pdf

# Count scatter plots
scatter_files <- list.files(results_path, pattern = "SuppFig7_scatterplot_.*\\.pdf", full.names = FALSE)
cat("✓ Scatter plots generated:", length(scatter_files), "files\n")

## ✓ Scatter plots generated: 16 files

cat("\nAll results saved to:", results_path, "\n")

## 
## All results saved to: /bmbl_data/yuzhou/collaborative/Sizun_lab/INDEPTH/SGCC/SGWT_results/DLBCL_GeoMX/GeoMXNew_code_data_publish/Results

sessionInfo()

## R version 4.4.2 (2024-10-31)
## Platform: x86_64-redhat-linux-gnu
## Running under: Red Hat Enterprise Linux 8.10 (Ootpa)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblaso-r0.3.15.so;  LAPACK version 3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] edgeR_3.42.4                limma_3.56.2               
##  [3] ImpulseDE2_0.99.10          msigdbr_7.5.1              
##  [5] DT_0.29                     knitr_1.44                 
##  [7] ComplexHeatmap_2.16.0       circlize_0.4.15            
##  [9] pheatmap_1.0.12             gridExtra_2.3              
## [11] ggrepel_0.9.6               SingleCellExperiment_1.22.0
## [13] SummarizedExperiment_1.30.2 Biobase_2.60.0             
## [15] GenomicRanges_1.52.0        GenomeInfoDb_1.36.3        
## [17] IRanges_2.34.1              S4Vectors_0.44.0           
## [19] BiocGenerics_0.52.0         MatrixGenerics_1.12.3      
## [21] matrixStats_1.0.0           ggpubr_0.6.1               
## [23] qs_0.25.5                   progress_1.2.2             
## [25] igraph_2.1.4                deldir_1.0-9               
## [27] lubridate_1.9.2             forcats_1.0.0              
## [29] stringr_1.5.2               purrr_1.1.0                
## [31] tibble_3.3.0                ggplot2_4.0.0              
## [33] tidyverse_2.0.0             tidyr_1.3.1                
## [35] readr_2.1.4                 GSVA_2.0.7                 
## [37] Seurat_5.3.0                SeuratObject_5.0.2         
## [39] sp_2.2-0                    dplyr_1.1.4                
## [41] Ternary_2.3.5              
## 
## loaded via a namespace (and not attached):
##   [1] spatstat.sparse_3.1-0    bitops_1.0-7             doParallel_1.0.17       
##   [4] httr_1.4.7               RColorBrewer_1.1-3       tools_4.4.2             
##   [7] sctransform_0.4.2        backports_1.5.0          utf8_1.2.6              
##  [10] R6_2.6.1                 HDF5Array_1.34.0         mgcv_1.9-1              
##  [13] lazyeval_0.2.2           uwot_0.1.16              GetoptLong_1.0.5        
##  [16] rhdf5filters_1.18.1      withr_3.0.2              prettyunits_1.1.1       
##  [19] progressr_0.14.0         textshaping_0.3.6        cli_3.6.5               
##  [22] Cairo_1.6-1              spatstat.explore_3.5-2   fastDummies_1.7.5       
##  [25] labeling_0.4.3           sass_0.4.7               S7_0.2.0                
##  [28] spatstat.data_3.1-6      ggridges_0.5.4           pbapply_1.7-2           
##  [31] systemfonts_1.3.1        dichromat_2.0-0.1        parallelly_1.45.1       
##  [34] rstudioapi_0.15.0        RSQLite_2.3.1            generics_0.1.4          
##  [37] shape_1.4.6              RApiSerialize_0.1.2      crosstalk_1.2.0         
##  [40] vroom_1.6.3              ica_1.0-3                spatstat.random_3.4-1   
##  [43] car_3.1-3                Matrix_1.7-4             fansi_1.0.4             
##  [46] abind_1.4-8              lifecycle_1.0.4          yaml_2.3.7              
##  [49] carData_3.0-5            rhdf5_2.50.2             SparseArray_1.6.2       
##  [52] Rtsne_0.16               blob_1.2.4               promises_1.2.1          
##  [55] PlotTools_0.3.1          crayon_1.5.2             miniUI_0.1.1.1          
##  [58] lattice_0.22-6           beachmat_2.16.0          cowplot_1.2.0           
##  [61] annotate_1.78.0          KEGGREST_1.40.0          magick_2.8.7            
##  [64] pillar_1.9.0             rjson_0.2.21             future.apply_1.11.0     
##  [67] codetools_0.2-20         glue_1.8.0               spatstat.univar_3.1-4   
##  [70] data.table_1.14.8        vctrs_0.6.5              png_0.1-8               
##  [73] spam_2.11-0              gtable_0.3.6             cachem_1.0.8            
##  [76] xfun_0.52                S4Arrays_1.6.0           mime_0.12               
##  [79] survival_3.7-0           RcppHungarian_0.3        iterators_1.0.14        
##  [82] ellipsis_0.3.2           fitdistrplus_1.1-11      ROCR_1.0-11             
##  [85] nlme_3.1-166             bit64_4.0.5              RcppAnnoy_0.0.21        
##  [88] bslib_0.5.1              irlba_2.3.5.1            KernSmooth_2.23-24      
##  [91] colorspace_2.1-0         DBI_1.1.3                DESeq2_1.40.2           
##  [94] tidyselect_1.2.1         bit_4.0.5                compiler_4.4.2          
##  [97] graph_1.78.0             DelayedArray_0.32.0      plotly_4.10.2           
## [100] stringfish_0.15.8        scales_1.4.0             lmtest_0.9-40           
## [103] SpatialExperiment_1.16.0 digest_0.6.33            goftest_1.2-3           
## [106] spatstat.utils_3.1-5     rmarkdown_2.24           XVector_0.40.0          
## [109] htmltools_0.5.6          pkgconfig_2.0.3          sparseMatrixStats_1.12.2
## [112] fastmap_1.1.1            rlang_1.1.6              GlobalOptions_0.1.2     
## [115] htmlwidgets_1.6.2        shiny_1.7.5              farver_2.1.2            
## [118] jquerylib_0.1.4          zoo_1.8-12               jsonlite_1.8.7          
## [121] BiocParallel_1.34.2      BiocSingular_1.16.0      RCurl_1.98-1.12         
## [124] magrittr_2.0.4           Formula_1.2-5            GenomeInfoDbData_1.2.10 
## [127] dotCall64_1.2            patchwork_1.3.2          Rhdf5lib_1.28.0         
## [130] Rcpp_1.1.0               babelgene_22.9           reticulate_1.32.0       
## [133] stringi_1.8.7            zlibbioc_1.46.0          MASS_7.3-61             
## [136] plyr_1.8.9               parallel_4.4.2           listenv_0.9.0           
## [139] Biostrings_2.74.1        splines_4.4.2            tensor_1.5              
## [142] hms_1.1.3                locfit_1.5-9.8           spatstat.geom_3.5-0     
## [145] ggsignif_0.6.4           RcppHNSW_0.4.1           reshape2_1.4.4          
## [148] ScaledMatrix_1.8.1       XML_3.99-0.14            evaluate_0.21           
## [151] RcppParallel_5.1.7       foreach_1.5.2            tzdb_0.4.0              
## [154] httpuv_1.6.11            RANN_2.6.2               polyclip_1.10-7         
## [157] clue_0.3-64              future_1.67.0            scattermore_1.2         
## [160] rsvd_1.0.5               broom_1.0.9              xtable_1.8-4            
## [163] RSpectra_0.16-2          rstatix_0.7.2            later_1.3.1             
## [166] ragg_1.5.0               viridisLite_0.4.2        memoise_2.0.1           
## [169] AnnotationDbi_1.62.2     cluster_2.1.6            timechange_0.2.0        
## [172] globals_0.18.0           GSEABase_1.68.0

---

Analysis completed successfully!

This comprehensive analysis generated all figures for Figure 5 and Supplementary Figure 7, including: - Pathway enrichment heatmaps for CD4T cells and macrophages - Ternary plots showing spatial correlations between cell types - Correlation analysis between immune dysfunction markers and EBV status - Interactive visualizations and statistical summaries